C $TITLE: 'FACTR'
C $NOFLOATCALLS
C
C
      SUBROUTINE FACTR(A,D,N,NDIM,IP,LD2)
C
C     SUBROUTINE TO FACTOR A MATRIX INTO A UNIT LOWER TRIANGULAR MATRIX
C     AND AN UPPER TRIANGULAR MATRIX USING THE GAUSS-DOOLITTLE ALGORITHM
C     PRESENTED ON PAGES 411-416 OF A. RALSTON--A FIRST COURSE IN
C     NUMERICAL ANALYSIS.  COMMENTS BELOW REFER TO COMMENTS IN RALSTONS
C     TEXT.    (MATRIX TRANSPOSED.
C
      INTEGER*4 N,NDIM
      INTEGER R,RM1,RP1,PJ,PR
      REAL*8 ELMAG,DMAX
CLARGE: A
      COMPLEX A
      COMPLEX*16 D,ARJ,SUM
      DIMENSION D(LD2)
C**
C**      DIMENSION IP(NDIM),A(NDIM,NDIM)
      DIMENSION IP(NDIM),A(NDIM,1)
C**
C     D      WRITE(*,*) '   FACTR: START N=',N,' NDIM=',NDIM
C**
      IFLG=0
C**
      SUM=DCMPLX(0.,0.)
      DO 9 R=1,N
C
C     STEP 1
C
      DO 1 K=1,N
      D(K)=A(R,K)
1     CONTINUE
C
C     STEPS 2 AND 3
C
      RM1=R-1
      IF (RM1.LT.1) GO TO 4
      DO 3 J=1,RM1
      PJ=IP(J)
      ARJ=D(PJ)
      A(R,J)=ARJ
      D(PJ)=D(J)
      JP1=J+1
      DO 2 I=JP1,N
      D(I)=D(I)-A(J,I)*ARJ
2     CONTINUE
3     CONTINUE
4     CONTINUE
C
C     STEP 4
C
      DMAX=DREAL(D(R)*DCONJG(D(R)))
      IP(R)=R
      RP1=R+1
      IF (RP1.GT.N) GO TO 6
      DO 5 I=RP1,N
      ELMAG=DREAL(D(I)*DCONJG(D(I)))
      IF (ELMAG.LT.DMAX) GO TO 5
      DMAX=ELMAG
      IP(R)=I
5     CONTINUE
6     CONTINUE
      IF (DMAX.LT.1.D-10) IFLG=1
      PR=IP(R)
      A(R,R)=D(PR)
      D(PR)=D(R)
C**
      SUM=SUM+A(R,R)
C**
C
C     STEP 5
C
      IF (RP1.GT.N) GO TO 8
      ARJ=1.D0/A(R,R)
      DO 7 I=RP1,N
      A(R,I)=D(I)*ARJ
7     CONTINUE
8     CONTINUE
      IF (IFLG.EQ.0) GO TO 9
      WRITE(*,10)  R,DMAX
      IFLG=0
9     CONTINUE
      RETURN
C
10    FORMAT (' PIVOT(',I3,')=',1P,D16.8)
      END
C
C
C
      SUBROUTINE SOLVES(A,B,Y,NEQ,NP,N,MP,M,IP,NRH,IFL1,IFL2,LD2,
     1 IRESRV)
C
C     SUBROUTINE SOLVES, FOR SYMMETRIC STRUCTURES, HANDLES THE
C     TRANSFORMATION OF THE RIGHT HAND SIDE VECTOR AND SOLUTION OF THE
C     MATRIX EQ.
C
      COMMON/SMAT/ SSX(16,16)
      INTEGER*4 NPEQ,NROW,NEQ,NP,N,MP,M
      INTEGER*4 IMAT,NPBLK,NLAST,NLSYM,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
      REAL*8 FNORM
CLARGE: A,B
      COMPLEX A,B
      COMPLEX*16 Y,SUM
      COMPLEX*16 SSX
      DIMENSION A(1),IP(1),B(NEQ,NRH),Y(LD2)
      COMMON/MATPAR/ICASE,NBLOKS,NPBLK,NLAST,NBLSYM,NPSYM,NLSYM,IMAT,
     1 ICASX,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
C**
C     D      WRITE(*,*) '   SOLVES: START IFL1=',IFL1,' IFL2=',IFL2
      IRESRV=IRESRV
C**
      NPEQ=NP+2*MP
      NOP=NEQ/NPEQ
      FNOP=NOP
      FNORM=1.D0/FNOP
      NROW=NEQ
      IF (ICASE.GT.3) NROW=NPEQ
      IF (NOP.EQ.1) GO TO 11
      DO 10 IC=1,NRH
      IF (N.EQ.0.OR.M.EQ.0) GO TO 6
      DO 1 I=1,NEQ
1     Y(I)=B(I,IC)
      KK=2*MP
      IA=NP
      IB=N
      J=NP
      DO 5 K=1,NOP
      IF (K.EQ.1) GO TO 3
      DO 2 I=1,NP
      IA=IA+1
      J=J+1
2     B(J,IC)=Y(IA)
      IF (K.EQ.NOP) GO TO 5
3     DO 4 I=1,KK
      IB=IB+1
      J=J+1
4     B(J,IC)=Y(IB)
5     CONTINUE
C
C     TRANSFORM MATRIX EQ. RHS VECTOR ACCORDING TO SYMMETRY MODES
C
6     DO 10 I=1,NPEQ
      DO 7 K=1,NOP
      IA=I+(K-1)*NPEQ
7     Y(K)=B(IA,IC)
      SUM=Y(1)
      DO 8 K=2,NOP
8     SUM=SUM+Y(K)
      B(I,IC)=SUM*FNORM
      DO 10 K=2,NOP
      IA=I+(K-1)*NPEQ
      SUM=Y(1)
      DO 9 J=2,NOP
9     SUM=SUM+Y(J)*DCONJG(SSX(K,J))
10    B(IA,IC)=SUM*FNORM
11    IF (ICASE.LT.3) GO TO 12
      REWIND IFL1
      REWIND IFL2
C
C     SOLVE EACH MODE EQUATION
C
12    DO 16 KK=1,NOP
      IA=(KK-1)*NPEQ+1
      IB=IA
      IF (ICASE.NE.4) GO TO 13
      I=NPEQ*NPEQ
      READ (IFL1) (A(J),J=1,I)
      IB=1
13    IF (ICASE.EQ.3.OR.ICASE.EQ.5) GO TO 15
      DO 14 IC=1,NRH
14    CALL SOLVE(A(IB),B(IA,IC),Y,NPEQ,NROW,IP(IA),LD2)
      GO TO 16
15    CALL LTSOLV(A,B(IA,1),Y,IP(IA),NPEQ,NEQ,NRH,IFL1,IFL2,LD2)
16    CONTINUE
C**
C     D     IF (NOP.EQ.1) WRITE(*,*) '   SOLVES: EARLY RETURN'
C**
      IF (NOP.EQ.1) RETURN
C
C     INVERSE TRANSFORM THE MODE SOLUTIONS
C
      DO 26 IC=1,NRH
      DO 20 I=1,NPEQ
      DO 17 K=1,NOP
      IA=I+(K-1)*NPEQ
17    Y(K)=B(IA,IC)
      SUM=Y(1)
      DO 18 K=2,NOP
18    SUM=SUM+Y(K)
      B(I,IC)=SUM
      DO 20 K=2,NOP
      IA=I+(K-1)*NPEQ
      SUM=Y(1)
      DO 19 J=2,NOP
19    SUM=SUM+Y(J)*SSX(K,J)
20    B(IA,IC)=SUM
      IF (N.EQ.0.OR.M.EQ.0) GO TO 26
      DO 21 I=1,NEQ
21    Y(I)=B(I,IC)
      KK=2*MP
      IA=NP
      IB=N
      J=NP
      DO 25 K=1,NOP
      IF (K.EQ.1) GO TO 23
      DO 22 I=1,NP
      IA=IA+1
      J=J+1
22    B(IA,IC)=Y(J)
      IF (K.EQ.NOP) GO TO 25
23    DO 24 I=1,KK
      IB=IB+1
      J=J+1
24    B(IB,IC)=Y(J)
25    CONTINUE
26    CONTINUE
C**
C     D      WRITE(*,*) '   SOLVES: RETURN AT END'
C**
      RETURN
      END
C
C
C
      SUBROUTINE SOLVE(A,B,Y,N,NDIM,IP,LD2)
C
C     SUBROUTINE TO SOLVE THE MATRIX EQUATION LU*X=B WHERE L IS A UNIT
C     LOWER TRIANGULAR MATRIX AND U IS AN UPPER TRIANGULAR MATRIX BOTH
C     OF WHICH ARE STORED IN A.  THE RHS VECTOR B IS INPUT AND THE
C     SOLUTION IS RETURNED THROUGH VECTOR B.    (MATRIX TRANSPOSED.
C
      INTEGER*4 N,NDIM
      INTEGER PI
CLARGE: A,B
      COMPLEX A,B
      COMPLEX*16 Y,SUM
      DIMENSION A(NDIM,1),B(NDIM),IP(NDIM),Y(LD2)
C**
C     D      WRITE(*,*) '   SOLVE: START'
C**
C
C     FORWARD SUBSTITUTION
C
      DO 3 I=1,N
      PI=IP(I)
      Y(I)=B(PI)
      B(PI)=B(I)
      IP1=I+1
      IF (IP1.GT.N) GO TO 2
      DO 1 J=IP1,N
      B(J)=B(J)-A(I,J)*Y(I)
1     CONTINUE
2     CONTINUE
3     CONTINUE
C
C     BACKWARD SUBSTITUTION
C
      DO 6 K=1,N
      I=N-K+1
      SUM=DCMPLX(0.,0.)
      IP1=I+1
      IF (IP1.GT.N) GO TO 5
      DO 4 J=IP1,N
      SUM=SUM+A(J,I)*B(J)
4     CONTINUE
5     CONTINUE
      B(I)=(Y(I)-SUM)/A(I,I)
6     CONTINUE
C**
C     D      WRITE(*,*) '   SOLVE: RETURN'
C**
      RETURN
      END
C
      SUBROUTINE LTSOLV (A,B,Y,IX,NROW,NEQ,NRH,IFL1,IFL2,LD2)
C
C     LTSOLV SOLVES THE MATRIX EQ. Y(R)*LU(T)=B(R) WHERE (R) DENOTES ROW
C     VECTOR AND LU(T) DENOTES THE LU DECOMPOSITION OF THE TRANSPOSE OF
C     THE ORIGINAL COEFFICIENT MATRIX.  THE LU(T) DECOMPOSITION IS
C     STORED ON TAPE 5 IN BLOCKS IN ASCENDING ORDER AND ON FILE 3 IN
C     BLOCKS OF DESCENDING ORDER.
C
      INTEGER*4 I2,I,J,K,NROW,NEQ,IDM1
      INTEGER*4 IMAT,NPBLK,NLAST,NLSYM,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
CLARGE A,B
      COMPLEX A,B
      COMPLEX*16 Y,SUM
      COMMON/MATPAR/ICASE,NBLOKS,NPBLK,NLAST,NBLSYM,NPSYM,NLSYM,IMAT,
     1 ICASX,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
      DIMENSION A(NROW,1),B(NEQ,NRH),Y(LD2),IX(NEQ)
C
C     FORWARD SUBSTITUTION
C
C     D      WRITE(*,*) '    LTSOLV: START'
      IDM1=1
      I2=2*NPSYM*NROW
      DO 4 IXBLK1=1,NBLSYM
      CALL BLCKIN (A,IDM1,I2,1,121,IFL1)
      K2=NPSYM
      IF (IXBLK1.EQ.NBLSYM) K2=NLSYM
      JST=(IXBLK1-1)*NPSYM
      DO 4 IC=1,NRH
      J=JST
      DO 3 K=1,K2
      JM1=J
      J=J+1
      SUM=(0.,0.)
      IF (JM1.LT.1) GO TO 2
      DO 1 I=1,JM1
1     SUM=SUM+A(I,K)*B(I,IC)
2     B(J,IC)=(B(J,IC)-SUM)/A(J,K)
3     CONTINUE
4     CONTINUE
C
C     BACKWARD SUBSTITUTION
C
      JST=NROW+1
      DO 8 IXBLK1=1,NBLSYM
      CALL BLCKIN (A,IDM1,I2,1,122,IFL2)
      K2=NPSYM
      IF (IXBLK1.EQ.1) K2=NLSYM
      DO 7 IC=1,NRH
      KP=K2+1
      J=JST
      DO 6 K=1,K2
      KP=KP-1
      JP1=J
      J=J-1
      SUM=DCMPLX(0.,0.)
      IF (NROW.LT.JP1) GO TO 6
      DO 5 I=JP1,NROW
5     SUM=SUM+A(I,KP)*B(I,IC)
      B(J,IC)=B(J,IC)-SUM
6     CONTINUE
7     CONTINUE
8     JST=JST-K2
C
C     UNSCRAMBLE SOLUTION
C
      DO 10 IC=1,NRH
      DO 9 I=1,NROW
      IXI=IX(I)
9     Y(IXI)=B(I,IC)
      DO 10 I=1,NROW
10    B(I,IC)=Y(I)
C**
C     D      WRITE(*,*) '    LTSOLV: RETURN'
C**
      RETURN
      END
